W : an alternative phenomenological coupling parameter for model systems 
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We introduce a parameter W(/3, L) = (it {\m\} 2 / {m 2 } — 2)/(tt — 2) which like the kurtosis (Binder 
cumulant) is a phenomenological coupling characteristic of the shape of the distribution p(m) of 
the order parameter m. To demonstrate the use of the parameter we analyze extensive numerical 
data obtained from density of states measurements on the canonical simple cubic spin- 1/2 Ising 
ferromagnet, for sizes L = 4 to L — 256. Using the W-parameter accurate estimates are obtained 
for the critical inverse temperature /3 C — 0.2216541(2), and for the thermal exponent v — 0.6308(4). 
In this system at least, corrections to finite size scaling are significantly weaker for the VK-parameter 
than for the Binder cumulant. 
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INTRODUCTION 

Studies of the critical properties of model systems us- 
ing numerical simulations are necessarily limited to sam- 
ples of finite size. Finite Size Scaling (FSS) techniques 
are essential in this context and a Renormalization Group 
Theory (RGT) of FSS is well established PK]. At the 
critical point the shape of the distribution of the order 
parameter p(m) (throughout we will use terminology ap- 
propriate to ferromagnetism) is independent of size L to 
within finite size correction factors. 

One widely used parameter characteristic of p(m) is 
the kurtosis of the distribution, U4 = (to 4 ) /(to 2 ) 2 . The 
kurtosis is often expressed in terms of the Binder cumu- 
lant |5] 

#I) = 1(3-^,I)) (1) 

because with this normalization g(0, L) = in the high 
temperature Gaussian limit (this is not strictly true for 
very small L) and g(po,L) = 1 in the low temperature 
ferromagnetic (non-degenerate ground state) limit (/3 = 
J/ksT is as usual the normalized inverse temperature). 
For a given sample geometry (such as a [hyper] cube) the 
thermodynamic (large L) limit of g(/3 c ,L) is a universal 
parameter for all systems in the same universality class. 
Again in the large L limit FSS theory shows that the 
slope dg{f3,L)/d(3 ~ L 1 ^ at (3 C where v is the standard 
thermal critical exponent. There are however corrections 
to FSS which must be taken into account at all finite L. 

The Binder parameter g(/3, L) is not the only distribu- 
tion " shape" parameter having these properties. We will 
introduce and illustrate on the simple cubic 5=1/2 Ising 
ferromagnet an alternative parameter W((3, L) which has 
some technical advantages at least in this case, in par- 
ticular having corrections to FSS which are weaker than 
those of the Binder parameter. If this turns out to be a 
general property the W-parameter could be very helpful 
for estimating critical properties numerically in more dif- 



ficult cases where simulations are intrinsically restricted 
to more moderate size samples. Already for the 3d Ising 
ferromagnet we obtain rather precise estimates for the 
critical inverse temperature /3 C and the thermal critical 
exponent v using this parameter. 

PHENOMENOLOGICAL COUPLINGS 

A "phenomenological coupling" is broadly a parame- 
ter which becomes L independent at f3 c in the thermo- 
dynamic limit [7j; as well as the Binder parameter the 
normalized second moment correlation length ^(/3.L)/L 
is an important phenomenological coupling. It should 
be noticed that below T c both the Binder parameter 
and £(/3,L)/L are denned in terms of non-connected 
distribution sums. While g((3, L) is defined such that 
g(oo,L) = 1, the conventional definition of £(/3, L) leads 
to £(00, L)/ L = 00 for a system with a non-degenerate 
ground state. A phenomenological coupling with a dif- 
ferent normalization, defined by R^((3, L) — £(/3, L)/(L + 
£(/3, L)), would be closer in spirit to the Binder cumulant. 

The standard RGT FSS expression with leading correc- 
tion terms for a phenomenological coupling R(/3, L) such 
as g(P,L), the normalized correlation length ^/L(/3, L), 
or the W(f3, L) to be introduced below, is 

R(f3, L) = R(u T L}l v ) + v u R^UrL 1 ^) L~ u + ■ ■ ■ 

« Rc + [dR/dr} c t tL x ^ + ■ ■ ■ + c u L~ u + ■ ■ ■ (2) 

where r is the thermal scaling variable (for instance 
t = 1 — /3/{} c ), u T is the thermal scaling field, cj = 8/v 
is a universal scaling correction exponent, and the other 
parameters (critical temperature and critical amplitudes) 
are non-universal constants appropriate for each partic- 
ular system. The second line is a good approximation as 
long as tL 1 ^ < 1. Thus 

R(j3 c ,L)=R c {l + c u L~" + ■■■), (3) 
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and 



[dR/d^^KL 1 '" (l + K u L- 



(4) 



Another important conclusion [8] is that the intersec- 
tion temperatures for R(/3,L) and R(j3,sL), denoted 
Across (L,s), converge as 



R — R ryj T,~( U+1 / V ) 

Pcross He 



(5) 



The parameter W(/3,L) which we introduce is a func- 
tion of the ratio of the variance of the modulus of m, 

Xmod = <(H ~ (H)) 2 ) = (to 2 ) - (|to|) 2 (6) 

to the variance of to 



X 



*» 2 ) 



(7) 



For a finite L ferromagnet in zero applied field, which is 
the case that we will discuss explicitly, the distribution 
p(m) is always symmetric so (m) = even below the 
critical temperature, thus x = (to 2 ) 

We will define the normalized parameter 



or 



W = 1 



w 



7T Xmod 
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where Ui = (|to|) 2 /(to 2 ). The normalization has been 
chosen such that, as for the Binder parameter, W(f3, L) — 

in the high temperature Gaussian limit and W(0, L) — 

1 in the low temperature ferromagnetic limit. As W(j3, L) 
is also a parameter characteristic of the shape of the dis- 
tribution p(m), it can be considered to be another "phe- 
nomenological coupling" and so will share all the formal 
finite size scaling properties of g(/3, L). 

By analogy with "kurtosis" (derived from the Greek 
word for a curve or bulge) we propose to name W the 
"dichokurtosis" , referring to the process of dividing into 
two parts, i.e. a unimodal distribution shifting into a 
bimodal one. 

As a demonstration, extensive data on W(f3, L) will be 
discussed for the canonical case of the simple cubic spin 
5 = 1/2 Ising ferromagnet. Though we will not discuss 
this point further here, the properties of the distribution 
p(|m|) are of particular interest when the regime T < 
T c is studied as well as T > T c . Above T c , (to) = 
in zero applied field; the connected and non-connected 
susceptibilities 



Xconn = (m ) - (m h ) 



and 



(10) 



(11) 



are identical. Below T c it is the connected susceptibil- 
ity which is physically significant in the thermodynamic 
limit. For finite L the distribution p(m) consists approx- 
imately of two peaks centered on ±(|to|). As (3 — j3 c and 
L increase this approximation gets better and better be- 
cause the peaks narrow so that 



<H> 



(m h ) 



(12) 



where (m/,) is the magnetization that would be measured 
in an infinitesimal applied field. Hence x, no( j(/3,L) ~ 

Xconn (A L). 

NUMERICAL METHODS 

The physical parameters for finite size samples from 
L = 4 up to L = 256 (16,777,216 spins) were esti- 
mated using a density of states function method. When 
studying a statistical mechanical model complete infor- 
mation can in principle be obtained through the density 
of states function. From complete knowledge of the den- 
sity of states one can immediately work with the micro- 
canonical (fixed energy) ensemble and of course also com- 
pute the partition function and through it have access to 
the canonical (fixed temperature) ensemble as well. The 
main problem here is that computing the exact density 
of states for systems of even modest size is a very hard 
numerical task. However, several sampling schemes have 
been given for obtaining approximate density of states, 
of which the best known are the Wang-Landau [5] and 
Wang-Swendsen [10] methods. In pT] various methods 
are discussed in a common framework. For work in the 
microcanonical ensemble the sampling methods give all 
the information needed. Using them one can find the 
density of states in an energy interval around the critical 
region and that is all that is required for most investi- 
gations of the critical properties of the model. It should 
be noted that there is no standard technique for estimat- 
ing the error bars in the outputs of this class of method, 
other than repeating the entire calculation a number of 
times which would be extremely laborious. 

For the present analysis a density of states function 
technique based upon the same method as in |12) was 
used though with considerable numerical improvements 
for all L studied here (adequate improvements to the 
L = 512 data set would unfortunately have been too 
time-consuming). The microcanonical (energy depen- 
dent) data were collected as described in [11]. We use 
standard Metropolis single spin-flip updates, sweeping 
through the lattice system in a type-writer order. Mea- 
surements take place when the expected number of spin- 
flips is at least the number of sites. For high tempera- 
tures this usually means two sweeps between measure- 
ments and four or five sweeps for the lower temperatures 
we used. Note that in the immediate vicinity of j3 c the 
spin- flip probability is very close to 50%. 
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For L = 256, the largest lattice studied here, we have 
now amassed between 500 and 3500 measurements on an 
interval of some 450000 energy levels, where most sam- 
plings are near the critical energy. For L = 128 we have 
between 5000 and 50000 measurements on some 150000 
energy levels. For L < 64 the number of samplings is of 
course vastly bigger. 

Our measurements at each individual energy level in- 
clude local energy statistics and magnetisation moments. 
The microcanonical data were then converted into canon- 
ical (temperature dependent) data according to the tech- 
nique in |13j . This gave us energy distributions from 
which we may obtain energy cumulants (e.g. the specific 
heat) and magnetization cumulants (e.g. the susceptibil- 

ity)- 

Typically around 200 different temperatures were cho- 
sen to compute these quantities, with a higher concen- 
tration near f3 c particularly for the larger L so that one 
may use standard interpolation techniques on the data 
to obtain intermediate temperatures. Magnetization dis- 
tributions p(m)(P,L) have also been obtained for sizes 
from L = 4 to L = 64. 
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FIG. 1: (Color online) Variation coefficient a/n for U2 (red 
circles) and U4 (blue squares) at /3 — 0.225 for L — 16 plotted 
versus the number of measurements n together with fitted 
lines with slope —1/2. The red line is 0.504/^/n and blue fine 
is 0.886/\Ai. The inset shows the ratio between these two for 
a range of j3, having the minimum 1.76 at /3 = 0.225. 



EQUILIBRATION TIMES 

We can make a critical comparison between the Binder 
parameter and the IT-parameter from the point of view 
of equilibration time. At the heart of the W-parameter 
is the ratio U2 — (\m\) 2 / (m 2 ) , just as the ratio U4 — 
(m 4 )/(m 2 ) 2 is the basis for the ^-parameter. Since the 
{/4-ratio involves a fourth moment we expect it to con- 
verge more slowly to its limit value than f/ 2 , which con- 
tains only a second moment. We have measured the 
speed of convergence by studying the respective varia- 
tion coefficients a/fi as a function of the number of mea- 
surements n. As usual a refers to the standard deviation 
of the measurements and \i to the average measurement. 
This allows us to compare the two, though the result 
will of course depend on the underlying distribution. We 
have chosen to look at this for a simple cubic lattice with 
L = 16 at ft = 0.225, a temperature slightly below where 
the distribution changes from unimodal to bimodal. 

We perform n measurements of |m|, m 2 and m 4 and 
take their respective averages, giving us estimates of 
(|m|), (m 2 ) and (m 4 ). The estimate of U2 is now sim- 
ply (\m\) 2 /(m 2 ) and for U4 we use (m 4 )/(m 2 ) 2 . Re- 
peating these n-estimates a number of times (75 times 
for n=100000 and 75000 times for n=100) gives us, in 
turn, an estimate of the variance a 2 of the U2- and U4- 
estimates. As \i we use the average U2- and t/4-estimates. 

In Figure [l] we show a j 1 \i versus n for U2 and U4 for 
the 3d-lattice with L = 16. We have fitted lines with 
slope —1/2 since we expect the variation coefficient to 
decrease at the rate l/y/n. We find that cr//i scales as 
roughly 0.504/^ for U 2 and 0.886/v^ for U A . 
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FIG. 2: (Color online) W((3, L) versus /3 for L = 4 (smallest 
slope) to L = 256 (strongest slope). The inset shows a zoomed 
in picture near j3 c . 



Squaring the factor 0.886/0.504 w 1.76 gives that U A 
requires 1.76 2 ~ 3.1 times as many measurements as U2 
to obtain the same statistical error a/fj, at (3 — 0.225. 

The factor 1.76 is actually close to a worst case scenario 
for this particular lattice. For higher temperatures, i.e. 
/3 < f3 c , this factor takes a value close to three and fol- 
lower temperatures, i.e. f3 > /3 C , the factor quickly ap- 
proaches a value close to four, Figure [l] inset. It also 
turns out that this worst case factor actually increases 
with L. For L — 8 we measured it to 1.74 at /3 = 0.23 
while for L = 32 we found it to be 1.83 at /3 = 0.2225. 
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FIG. 3: (Color online) dW/d/3 versus /3 for L = 
32, 64, 128, 256 where the maximum increases with L. The 
red vertical line is located at /3 C = 0.2216541. 
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FIG. 4: (Color online) Crossing points /3 C ross versus 1/i 2 ' 40 
for L = 6,8, 16, 32, 64, 128 for W (red circles) and g (blue 
squares) and fitted lines. 
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FIG. 5: (Color online) Log-log plot of dW/d/3 (red circles) 
and dg/d/3 (blue squares) versus L at f3 = f3 c = 0.2216541. 
Lines have slope 1/u. 
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FIG. 6: (Color online) L~ 1/v dW/d/3 - 0.32 (red circles) and 
L~ 1/v dg/dP (blue squares) at P = p c = 0.2216541 versus 
L. The W^-points have been translated by —0.32 for easier 
comparison. 
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FIG. 7: (Color online) W(/3 C ,L) + 0.22 (red circles) and 
g(p c ,L) (blue squares) versus 1/L U with uj — 0.814 for 
L = 6,8, 16, 32, 64, 128, 256. The W-points have been trans- 
lated by 0.22 for easier comparison. 



THE 3D ISING FERROMAGNET 

The spin- 1/2 Ising ferromagnet on a simple cubic lat- 
tice is an archetypical model system which has been very 
extensively studied. Although there exist no exact re- 
sults for any of the critical parameters, /3 C and the criti- 
cal exponents are known to high precision thanks to RGT 
theory, high temperature series expansions (HTSE), and 
numerical simulations (see Refs. [l4TH6j ). There is con- 
sensus that for this system j3 c s» 0.221655, and for the 
universality class v sa 0.630, and to « 0.81. We will test 
our W-data against these values shortly. 

We show in Figure [2] an overall view of the behavior of 
W(j3,L). It can be seen that on the scale of the figure 
the curves W(/3, L) appear to intersect at a unique L- 
independent inverse temperature which can obviously be 
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identified with f3 c . The derivatives dW/df3 peak strongly 
at PwmaAL) which at large L approaches /3 C , see Fig- 
ured 

A blow-up of W(/3, L) in the critical region, Figure [2] 
inset, shows that there are finite size corrections lead- 
ing to a weak size dependence of the [W(/3, L), W((3, 2L)] 
crossing points. A plot of the intersection temperatures 
Across versus l/L UJ+1 / lJ , where w + \jv w 2.40, is shown 
in Figure [4] for both W and g (the points for L = 64, 128 
are not visible due to the fast data collaps). Fitting a 
straight line to the W-points for L > 6 versus 1/L W+1 / I/ 
gives f3 c — 0.2216541(5). We have here excluded the 
point L = 4 since it appears to deviate from the others 
in this case. The error estimate is based on how the re- 
sult depends on excluding a point from the fit and on 
allowing the exponent uj + 1/v to take different values 
between 2.39 and 2.41. In fact, a best fit of the cross- 
ing points for L > 6 to a simple formula Cq + c x L~ x 
gives on average, taken over fits after excluding one point, 
c « 0.2216540(3) and A = 2.40(3), where the error esti- 
mates correspond to the standard deviation of the data 
set. Since v = 0.630 is known to a higher precision we 
therefore get to — 0.814(30) which agrees with previous 
estimates. 

Both of our /3 c -estimates are consistent with the most 
precise values from standard Monte-Carlo simulations 
P c = 0.22165452(8) [15], /? c = 0.2216546(3) [12], and 
from high temperature series analyses, f3 c — 0.221655(2) 
[15j . However, the g-data seem to require more correc- 
tion to scaling than the H'-data. If we want to fit a line 
to the crossing points for g versus L w+1 / V then we need 
to drop two more points (L — 6, 8) to get anything like 
this precision on a /3 c -estimate. 

Henceforth setting f3 c — 0.2216541, let us proceed to 
investigate the derivative data [dW/d(3)fj c and [dg/dj3)j3 c 
against L, which are shown as log-log plots in Figure [5] 
The slopes should be equal to 1/v in the large L limit. It 
can be seen that both series of points lie close to 1/0.63 
(slope of the lines). 

Let us make a more demanding analysis of the slopes 
1/v by fitting lines to fc-subsets of the points. Since we 
have 9 data points, i.e. we use L > 4, each k then gives 
us ( fe ) different slopes. If the data show any sign of in- 
consistency or a dependency on L then we expect this 
to show up in the form of different medians and/or dif- 
ferent slope intervals. However, we get v = 0.6308 for 
k = 3, . . . , 9, with the same value for both median and 
mean. The quartile deviation of each slope set is about 
0.0004 for k = 4, . . . , 7. We therefore receive the esti- 
mate v — 0.6308(4). It should be noted that only for the 
last three points of the <?-data do we receive a slope that 
agrees with this estimate. 

An alternative way of locating j3 c is to locate the tem- 
perature j3 c where the scaling of the derivatives depend 
least on different L. Choosing e.g. subsets of size k = 4 
the narrowest set of slopes is obtained for f3 c — 0.2216541, 



give or take a step or two in the last decimal. Since 
this agrees with our previous two estimates of /3 C we can 
now give our final estimate of the critical temperature as 
/3 C = 0.2216541(2). 

Having established /3 C and v we plot the derivatives of 
W and g in the more demanding form [dW/d^p^/L 1 ^ 
and [dg/d^pjL 1 ^ in Figure^ The g-data clearly show 
characteristic FSS corrections 

g(/3 c , L) = g{(5 c , oo) (l + a u + ) (13) 

at small and moderate L while the VK-data show only 
weak and apparently random scatter due to statistical er- 
rors, i.e. the analogous correction term for W((3 C , L) ap- 
pears negligible within the present precision. This means 
that to extract an estimate of v a two parameter fit is 
sufficient for the W derivative data while a four param- 
eter fit is needed for the g-data. This is important as 
it means that at least in the present case the estimates 
from W(f3 c ,L) are intrinsically more precise. 
It was estimated in [TB] that 

g(j3 c , L) = 0.69778(13) (l + 0.1788(36) £-°- 82 ( 3 ) + . . . ) 

(14) 

and our g-data, Figure [6] are in excellent agreement with 
this correction factor for g(f3 c ,L). We estimate the crit- 
ical values to be W(/? c ,oo) = 0.468(2) and <?(/3 c ,oo) = 
0.697(2), see Figure [7j where the error stems from which 
points are excluded from the fit. The value for g agrees 
with the formula above but the accuracy is not as good. 
Also, we would like to mention that at the temperature 
where the magnetisation distribution shifts from uni- 
modal to bimodal, i.e. where [d 2 p(m)/dm 2 ] Q — 0, we 
found the asymptotic value of W to be about 0.208 and 
for g it takes a value near 0.433. 

There are already many accurate estimates of v for 
the 3d Ising universality group. Renormalization group 
studies [g give v = 0.6304(13) and v = 0.6305(25). The 
main difficulty concerning either HTSE or MC analyses 
lies principally with the problem of properly allowing for 
corrections to scaling. The amplitudes of the corrections 
vary from system to system, favorizing meta-analyses of 
data on many systems in the same class. Butera and 
Comi [TS] obtain v = 0.6299(2) from a global analy- 
sis of HTSE data for Ising ferromagnets with spin S 
running from 1/2 to oo on both sc and bec lattices, 
all systems lying in the same universality class. Their 
sc S = 1/2 HTSE results standing alone were consis- 
tent with this value but were less accurate (0.632(2) 
or 0.6277(30) depending on the analysis method used). 
Deng and Blote TB] obtain an entirely independent global 
estimate v = 0.63020(12) from simultaneous Monte Carlo 
analyses on a set of eleven systems all in the same uni- 
versality class. It is gratifying that the present results 
on one single system are consistent with and practically 
as accurate as these global "best estimates" from HTSE 
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and MC. It would be interesting to establish whether 
the weak FSS correction for [dW(P,L)/d0\p e is a gen- 
eral property or is specific to this particular system. 

CONCLUSION 

We introduce an alternative distribution "shape" pa- 
rameter W(P, L) for numerical studies of the critical 
properties of model systems. As an illustration we use 
this parameter in an analysis of extensive data sets ob- 
tained through a density of states technique applied to 
simple cubic 5 = 1/2 Ising ferromagnet samples of size 
up to L — 256. In this system at least, corrections to scal- 
ing for W(j3 c ,L) are considerably weaker than those for 
the canonical Binder cumulant ,g(/3 c , L) and the equilibra- 
tion time to obtain data to a similar degree of precision is 
significantly lower. We obtain estimates for the critical 
inverse temperature f3 c — 0.2216541(2) and the critical 
exponents v = 0.6308(4) and oj = 0.814(30), based only 
on W-data, which are compatible with and almost as ac- 
curate as values from previous Monte Carlo |16) and high 
temperature series expansions |15j . 
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